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Abstract. We explore a mechanism of pattern formation arising in processes described 
by a system of a single reaction-diffusion equation coupled with ordinary differential 
equations. Such systems of equations arise from the modeling of interactions between 
cellular processes and diffusing growth factors. We focused on the model of early car- 
^•f-) cinogenesis proposed by Marciniak-Czochra and Kimmel, which is an example of a wider 

i-H class of pattern formation models with an autocatalytic non-diffusing component. We 

present a numerical study showing emergence of periodic and irregular spike patterns 
due to diffusion-driven instability. To control the accuracy of simulations, we develop 
a numerical code based on the finite element method and adaptive mesh. Simulations, 
t— h supplemented by numerical analysis, indicate a novel pattern formation phenomenon 

<^ based on the emergence of nonstationary structures tending asymptotically to the sum 

of Dirac deltas. 

Key words: diffusion-driven instability, spike patterns, numerical simulations, reaction- 
diffusion equations, mass concentration. 
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1. Introduction 

Classical mathematical models of biological or chemical pattern formation have been 

i i developed using reaction-diffusion equations, see eg. [3J EJ QUI EH E3] and references 

therein. In that framework there exist essentially two mechanisms of formation of stable 
j> spatially heterogeneous structures, 

• diffusion- driven instability (DDI) which leads to destabilization of a spatially ho- 
(*C) mogeneous steady state and emergence of Turing patterns, 

• a mechanism based on the multistability and hysteresis in the kinetic system which 
allows for the formation of transition layer patterns far from equilibrium. 

Both mechanisms can also coexist yielding a complex dynamics of the system as, for 
example, in the Lengyel- Epstein model of chemical reactions [HI 123] • 
£> The Turing phenomenon is related to a local behavior of solutions of a reaction-diffusion 

system in the neighborhood of a constant solution that is destabilized via diffusion. Pat- 
terns arising through a bifurcation can be spatially monotone or spatially periodic. The 
mechanism responsible for such behavior of model solutions is called a diffusion-driven 
instability (Turing- type instability), which can be formulated in the following way. 

Definition 1.1 (Diffusion-driven instability (DDI)). A system of reaction- diffusion equa- 
tions exhibits DDI (Turing instability) if and only if there exists a constant stationary 
solution which is stable to spatially homogeneous perturbations, but unstable to spatially 
heterogeneous perturbations. 

The original idea was presented by Turing on the example of two linear reaction- 
diffusion equations [21 J. Due to the local character of Turing instability, the notion has 
been extended in a natural way to the nonlinear equations using linearization around a 
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constant positive steady state. However, nonlinear systems may have multiple constant 
steady states yielding existence of heterogeneous structures far from the equilibrium. In 
such cases, global behavior of the solutions cannot be predicted by the properties of the 
linearized system and a variety of possible dynamics depending on the type of nonlin- 
earities can be observed. On the other hand, Turing instability can be exhibited also in 
degenerated systems such as reaction-diffusion-odes models or integro-differential equa- 
tions, for example shadow systems obtained through reduction of the reaction-diffusion 
model, [7], [IT]. Following all these observations and the character of Turing's original 
system, we define Turing patterns in the following way. 

Definition 1.2. By Turing patterns we call the solutions of reaction- diffusion equations 
that are 

• stable, 

• stationary, 

• continuous, 

• spatially heterogeneous and 

• arise due to the Turing instability (DDI) of a constant steady state. 

Recently, it has been shown that if DDI property is exhibited by a system of a single 
reaction-diffusion equation coupled to an ordinary differential equation with autocatal- 
ysis of non-diffusing component. Then, it does not lead to Turing patterns, namely all 
continuous patterns are unstable [IS]. As a consequence the question for the long-term 
behavior of solutions arises. It has been previously shown that a diffusion- driven blow-up 
in systems of react ion- diffusion equations can occur in finite time, [22] . Even more, blow- 
up in finite time in L°°, but global existence of weaker solutions has been shown, leading 
to so called "incomplete blow-up", see e.g. [18] for uniform boundedness in L 1 . 

In the current paper we present a phenomenon of diffusion- driven unbounded growth 
and formation of dynamic spike pattern converging asymptotically to a sum of Dirac 
deltas. For a reaction-diffusion-ode model arising from applications in biology, we show 
that introducing diffusion in the ODEs system not only destabilizes the constant steady 
state, but also leads to an unbounded growth of model solutions. Since the solutions of 
the system with zero diffusion are uniformly bounded, we call the observed phenomenon 
the diffusion- driven unbounded growth. The total mass (L 1 norm) of the solutions is 
uniformly bounded but it concentrates in isolated points for time tending to infinity. 
Using numerical simulations, we investigate how the shape of emerging patterns depends 
on initial conditions and the scaling coefficient (size of diffusion versus domain size). 
Interestingly, we find out that the shape of observed patterns are superposition of a 
near-equilibrium effect of diffusion-driven instability and a far-from-equilibrium effect of 
multistability exhibited by the model. 

2. Problem formulation 

We study a reaction-diffusion-odes model of the diffusion-regulated growth of cell pop- 
ulation, which has the form of two ordinary-partial differential equations 



(2.2) 



(2.1) 





for x e [0,1], t > 0, 
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supplemented with homogeneous Neumann (zero flux) boundary conditions for the func- 
tion w = w(x, t) 

(2.3) w x {0, t) = w x {l, t) = for all t > 0, 
and with nonnegative initial conditions 

(2.4) u(x, 0) = Uq(x), w(x, 0) = wq(x). 
ax, dx, D w , Kx denote positive constants. 

In this paper we focus on one-dimensional domain [0, 1] for a clarity of presentation. 
The results can be obtained also for a model defined on two-dimensional space domain. 
Obviously, in such case a structure of spatial patterns is richer. Nevertheless, the main 
aspect of the pattern formation phenomenon exhibited by this model, i.e. evolution of 
spike patterns of mass concentration, is preserved independent on the dimension of the 



spatial domain. Model (2.1)-(2.4) is a rescaled reduction of the model 



(2.5) u t = [a d c )u, for iG [0,1], t > 0, 

\ u + v J 

(2.6) v t = —d b v + au 2 w — dv, for x G [0,1], t > 0, 

(2.7) w t = —w xx — d g w — au 2 w + dv + k, for x 6 (0, 1), t > 0, 

supplemented with homogeneous Neumann (zero flux) boundary conditions for the func- 
tion w = w(x, t). 

Model ( 2.5[ )-( [2~7 ) was proposed in [IT] as a receptor-based model of spatially distributed 



growth of a clonal population of pre-cancerous cells and its extensions and modifications 
were studied in [T2], [13] . The reduction was proposed in [5], but without further numerical 
or analytical investigation. 

In case of the spatial domain being the unit square, approximation of solutions of model 



(2.5)-(2.7) have been performed in [5]. Numerical simulations of the models showed qual- 
itatively new patterns of behavior of solutions, including, in some cases, a strong depen- 
dence of the emerging pattern on initial conditions and quasi-stability followed by rapid 
growth of solutions. However, recently it has been shown using linear stability analysis of 
nonconstant steady states that all stationary solutions of this model, both continuous and 
discontinuous, are unstable [H]. A question arises if the model exhibits a formation of 
any pattern, which persist for long times. Our present research is focused on understating 
these phenomena and answering questions on pattern formation in such class of models. 

3. Analytical results 



In the remainder of this paper, we consider the system (2.1 )-(2.4). It has been obtained 
using a quasi-stationary approximation assuming that the dynamics of v variable is faster 
than the dynamics of other variables. In the present paper, we focus on the reduced model, 
since it is the simplest reaction-diffusion-ode model exhibiting the spike pattern formation 



mechanism. A rigorous link between the solutions of the original model (2.5)-(2.7) and 
its two-equations approximation has been recently shown in 



3.1. Existence of solutions. Existence of global, classical solutions can be proven within 
the framework of ordinary differential equations and the theory of linear semigroups, see 
e.g. [19, 20] . Moreover, it can be shown using maximum principle that the solutions 
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remain positive for positive initial conditions. 



3.2. Existence of steady states. The analytical results concerning existence of regular 



stationary patterns of (2.1 )-( 2.2 ) can be summarized in the following theorem: 



Theorem 3.1. Under assumptions a\ > d\ and K\ > 2 , system (2.1)-(2.2) has the 
following smooth stationary solutions 

• constant steady states (u Q ,w ) = (0, ki), (u + ,w+) = f + yj(f ) 2 - (^^Y 

and (u-,w~) = ( a ^ dl *^2~\J{^t) 2 ~ (^^d^) 2 ) being stationary solutions of the 
kinetic system. 

• a unique strictly increasing solution W and a unique strictly decreasing solution 
W; U is defined by U = ^ dl ^. 

• a periodic solution W with n modes, increasing on intervals [0, -] and its sym- 
metric counterpart W(x) = W n (l — x), where n £ N depends on the diffusion 
coefficient; and the periodic function W £ C([0, 1]) is defined in the following 

W(x-m dla xef^^l 

V n ) L n ' n J 

W (W_ X \ dla XE r22±l 32±2] 
for every j G {0, 1, 2, 3, ...} such that 2j + 2 < n. U is defined by U = ^ . 



The proof of this statement is deferred to the Appendix. 
3.3. Stability of steady states. We investigate stability of the solutions described in 



Theorem 3.1, item (i): 

The operator resulting from linearization of (2.1 )-( 2.2 ) around ( ~, w) reads in the 
matrix-form: 



(3.1) 



J :-- 



an 0,12 
«2i a 2 2 + D W A 



a\ — d\ 



ai 



di 



dl 



Assume that a solution of 4 



oi— 1 V (ai — di)ui / 

J0 with homogeneous Neumann boundary conditions is 
of the form = where 0^ denotes the eigenvector of the Laplace operator associated 
to the kth eigenvalue. Then, the dispersion relation, i.e. the dependence of eigenvalues 
of the problem linearized at a constant steady states with the eigenvalues of the Laplace 
operator, see Fig. 3.1 , is defined by 

(3.2) disp(A,A;) = det 



an — A 

0-21 a 22 



Ol2 

D w (irk) 2 - X 



where (nk) 2 is the k-th eigenvalue of the Laplace operator considered on C 2 (0, 1). There- 
fore, A is an element of the point spectrum of J if disp(A, k) — for A ^ an, 0. 

Proposition 3.2. Under assumptions ai > di and n\ > 2 a ^ di , the following holds 

• {uq,wq) is a stable stationary solution of (2.1)-(2.2) and its kinetic system. 

• (u + ,w + ) is an unstable stationary solution of (2.1)-(2.2) and its kinetic system. 

• is an unstable stationary solution of (2.1)-(2.2). 
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• (tt_,w_) is a stable stationary solution of the kinetic system of (2.1)-(2.2) if and 
only if at least one of the following conditions is satisfied: 

1) k\>2 
/ 1 ai(ai— dj) 

-2 J Ol > j — K ? > — 77 77. 

/ x ai — 1 J. oi oi— oi(oi— aij 

Additionally, there exist infinitely many positive eigenvalues of the operator resulting from 



a linearization of (2.1)-(2.2) at (u-,W-) and ai ai d\ and — oo are their only limit points. 



The proof of this proposition can be found in the Appendix. 



A solution with initial conditions close to (u 0} w ) is shown in Figure 7.4 





Figure 3.1. Roots of the dispersion relation for the operator resulting from 
a linearization of (2.1)-(2.2) around The parameters are a\ = 2,di = 

1, k% = 3, D w = 2. left: A + . right: A_. We see that there exist infinitely many 



positive eigenvalues. 

Moreover, all steady states except (u , w 
from Theorem 2.1 and Corollary 2.7 in 



are linearly unstable. The latter results directly 
due to autocatalysis of u for all u, w > 0. 



Theorem 3.3. Under assumptions a\ > d\ and K\ > 2 a ^ d , all steady states (U,W) G 
L°°(0, 1) x H l (0,l) of (2.1)-(2.2) with U > on a non-zero-measure set are linearly 
unstable. 



3.4. Boundedness properties. Numerical solutions of system (2.1)-(2.2) presented in 



the following chapter show unbounded growth of spikes. To understand the underlying 
phenomenon, we summarize here results on the boundedness of solutions of the model 
with and without diffusion. It can be easily shown that the mass of solutions of system 



(2.1 )-(2.2) is uniformly bounded in time, see Lemma 3.4 Therefore, unbounded growth of 



the solutions may happen at most in isolated points of the spatial domain. Furthermore, 
to exclude blow-up induced by unbounded solutions of the kinetic system, such as shown 



in [TJ, we check boundedness properties of the kinetic system, see Lemma 3.5 
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Lemma 3.4. Let u(x,t),w(x,t) denote a solution of (2.1)-(2.2) for positive initial con- 
ditions. Then, it holds 



(3.3) 

(3.4) 
(3.5) 



1 K\ 

limsup (— \\u(t)\\ L1 + \\w{t)\\ L1 ) < 
t->oo Oj\ mm(ai, 1, 

ax 



limsup ||u(*)|| L i < 
t^oo mm(ai, l, 

limsup || w(t) \\ L1 < K\. 



Ku 



Moreover, the solution of the kinetic system of ( 2.1[ )-(2.2) is uniformly bounded in time 



Lemma 3.5. Let u(x,t),w(x,t) denote a solution of the kinetic system of (2.1) -(2. 2) for 
positive initial conditions. Then holds 



(3.6) 
(3.7) 



limsup u(t) < 
t^oo mm(di,l; 

limsupw(t) < K\. 

t—>oo 



Both lemmas can be proven similarly as it was shown in [H] for the three equation 
model. More details are deferred to the Appendix. 

We conclude from Lemma [3.51 that the mass concentration observed in numerical sim- 
ulations does not result from a blow-up of the solution of the kinetic system. 

4. Numerical approach 



Numerical approximations of solutions to system (2.1 )-( 2.2 ) presented in this paper are 
obtained using the program library deal.ii, [2J. 

Simulations using adaptive grid refinement based on cell-wise evaluation of the proposed 



error indicators in |3j show a growth of spikes, see Figure 5.2 



The question of what is seen in numerical simulations motivated us to undertake a numer- 
ical study of the pattern formation phenomenon. To allow a rigorous argumentation using 
classical finite-element analysis, we investigate the asymptotic behavior of the numerical 
solution for spatially homogeneous meshes. 

For space discretization, we use a finite-element scheme with piecewise linear, globally 
continuous ansatzfunctions. The time discretization is performed using the implicit Euler 
scheme or the Crank-Nicholson scheme. 



Convergence for such scheme for solutions of systems of type (2.1)-(2.2) is well known, 
see [3]. 

The low order of the space discretization is due to the fact that preliminary simulations 
already showed emergence of spikes, corresponding to a large second derivative in space. 

5. Numerical analysis of the pattern formation phenomenon 
We choose parameters 
(5.1) a± = 2, d\ = 1, K\ = 3 

and diffusion coefficient D w = 6. 

A numerically obtained solution for different parameters, a\ = 2.5, d\ = 1.5, K\ = 4 is 



shown in the Appendix in figure 7.3 It shows qualitatively the same behavior. We recall 



that system (2.1)-(2.2) exhibits Turing type diffusion driven instability, but all positive 



steady states are linearly unstable. 
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5.1. Unbounded growth and spike formation. Initial conditions are chosen as a 



perturbation of the stable stationary solution (u,w) of the kinetics system of (2.1)-(2.2): 



(5.2) 



U (x) = U_ + €ip(x) 
Wq(x) = W-, 



where the 'perturbation function' p(x) satisfies the following conditions: 

(5.3) p is a polynomial of degree two on (0, s — e), (s — e, s + e), (s + e, 1) 

(5.4) peC\0,l), 

(5.5) p'(0)=p'(l) = 0, 

(5.6) p(0)=p(l) = -l, 

(5.7) p(s) = 1, 



and is thereby uniquely defined by the pair (s, e). The explicit formula for p can be 



found in the appendix, (7.33), an illustration can be found in figure 5.1 In figure 5.2, the 




0.2 0.4 0.6 0. 

x 



Figure 5.1. Illustration of the perturbation function p(x), defined by 



(5.3)-(5.7) for s = 0.4, e = 0.1. max x& np(x) is always assumed in 
(s - e, s + e). 



solution for initial conditions (5.2) for s = 0.4, ei = 0.05, e = 0.1 is shown. We observe 



exponential growth in a single point and decay towards zero otherwise. The maximum 
value of the numerically obtained solution keeps growing. 



Spike position and initial conditions. Simulations performed using the param- 
2, d\ — 1, K\ — 3, D w = 6 and initial conditions (5.2) show emergence of spikes 



5.2 

eters d\ 

at the maximum of the initial conditions, see Table [TJ For the linearized problem, this is 
heuristically reasonable since almost all eigenmodes of the Laplace Operator are unstable 



with almost the same eigenvalue, see Lemma 7.2 or figure 34] for an illustration. 
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Figure 5.2. Numerically obtained solution for initial conditions (5.2) with 
s = 0.4, e = 0.1, ei = 0.05 and diffusion coefficient D w — 6, u = 

2.62, w = 3 ~ 2 V ^ ~ 0.382. Left: component u. Right: component w. We 
observe formation of a spike at x ~ 0.43, which keeps growing exponentially 
in time. 



Lemma 5.1. Let J denote the operator resulting from the linearization of (2.1)-(2.2) 
around («_,«;_) and consider the initial value problem 



(5.8) 



d 
db 



1> 



J 



with homogeneous Neumann boundary (zero-flux) conditions for ip. 

As initial conditions take (ipo, Po) = (<f>k, \ - a2 2+D ^k 2 ^' w here <fik denotes the eigen- 
function of the Laplace- Operator (Neumann) associated to the kth eigenvalue and (o^) 
denotes the Jacobian of the kinetics system at wJ). 

Then e A+ ^^(^o,Po) is the solution of (5.8) for homogeneous Neumann boundary condi- 
tions. 

Heuristically speaking, the initial perturbation of u is self-amplifying for large D w if it 
is much larger than the perturbation of w, because X+(k) —> a u > 0. 
This heuristic implication leads to the question what happens for more complex initial 
conditions. In Figure |5.3[ we plot the numerical solution for initial conditions 



(5.9) 



Uq(x) = U- — ecos(47nr) 

Wq(x) = 



and in figure 5.4 the numerical solution for initial conditions 

Unix) = U- — e cos(47nr 2 ), 
(5.10) _ K h 

w {x) = w_. 
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Figure 5.3 
diffusion coefficient D 



Numerical solution for initial conditions (5.9) with e 

2 



it 



2.62, w 



0.05 and 
0.382. Left: 



component u. Right: component w. We observe formation of spikes at the 
position of local maxima of the initial conditions. 




Figure 5.4. Numerical solution for initial conditions (5.10) with e 



and diffusion coefficient D,, 



2, u 



3-V5 



0.05 
0.382. Left: 



- 2.62. w , 

component u. Right: component w. We observe formation of spikes at the 
position of local maxima of the initial conditions and faster growth at x = ^ 
than at x = |. 

We note that the initial perturbation seems to be indeed self amplifying. This does 
not explain the long-time behavior, but numerical simulations indicate that spikes grow 
close to the maxima of the initial conditions. We also note for initial conditions (5.10) 
that the spike for larger x grows faster. Real and imaginary part of the numerically ob- 
tained Finite Fourier Transform f(co) := J* cos(4irx' 2 )e t7ruJx dx and the growth rate of the 
perturbation at the maxima of uq are shown in Figure 7.1 , Figure 1_J2_ shows the growth 
rate of perturbation (5.9). 
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5.3. Varying the diffusion coefficient. The roots of the dispersion relation have the 
form 



(5.11) \±(k 2 



tr(A) - (irk) 2 D, 



± 



tr(A) - {-Kk) 2 D % 



\A\ + {irk) 2 D w a u 



We know that X-(k ) — > — oo and X + (k ) — > an > as k — > oo, see Lemma 7.1, 7.2 



Additionally, it holds A + (0) < since is a stable steady state of the kinetic 

system of (|2TI])-([272]) . 

It follows that there exist stable eigenmodes of the Laplace Operator, because X-(k 2 ) < 
and 



(5.12) 



\+{k 2 ) < k 2 < 



a u rr 2 D w 

This implies dampening of the low frequency part of the initial perturbation (4>,tp). 
First, we choose the same initial conditions and parameters as in figure 5.2, but vary the 
diffusion coefficient D,,,. 




3+V5 
2 



Figure 5.5. Numerical solution for initial conditions (5.2) with s 
0.1,62 = 0.05 and diffusion coefficient D w — 1, u 
0.382. Left: component u. Right: component w. 



2.62, w 



0.4, e 



For smaller diffusion coefficient 
initial conditions, see Figure 



5.5 



D w = 1, we observe growth of multiple spikes for the same 
We observe a self-amplification of the high-frequency 
part of the initial perturbation. The short-time behavior is therefore similar to the idea 
of a 'dominant' eigenvalue in classical Turing type models. However, in this case, we can 
speak of a 'self amplification of the part of the initial perturbation with sufficiently high 
wavenumber'. 
We define, 
(5.13) 



D 



w,k 



1 -Ad\ + (ai - di) 2 K 2 + k(cli - di)^jK\(ai - di) 2 - Ad{ 



an(-Kk) 2 (nky 



2d 2 



Figure 5.6 shows the solution for the corresponding D W) \ for model (2.1 )-(2.2). Table [2] 



shows the number of spikes for further variation of 
s = 0.4, e = 0.1, e 1 = 0.05. 



for initial conditions (5.2) with 
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Figure 5.6. Numerical solution for initial conditions (5.2) with s 
0.1, ei = 0.05 and diffusion coefficient D w = 5.8541 »s -D TO ,i, u - 
2.62, w = 3 ~ 2 V ^ w 0.382. Left: component w. Right: component u>. 



0.4, e 

3+v^ 



5.4. Evolution of mass. Our simulations indicate a growth of one or multiple spikes, 
u(x) — > oo for some x as t — > oo, and decay in all other x. We therefore investigate 



the evolution of the L 1 -norm of u. Figure 5/7 shows the evolution of mass of the solu- 
tion shown in figure 52 for homogeneous spatial mesh size h = 2 -16 and homogeneous 
temporal mesh size k = 2.5 ■ 10~ 4 . The convergence order is shown in the appendix, 
see Figure 7.5 7.6 Lemma 3.4 states that the the mass of the solution, \\u( 



Lemma 3.4 states that the the mass of the solution, ||w(£)|| L i is uni- 
formly bounded. However, an important question when modeling natural phenomena is 
positivity of mass if there is no extinction. The numerical simulations of the evolution 
of the mass suggests that it stays strictly positive. Therefore, based on the numerical 
simulations we have a conjecture that the solutions converge asymptotically to the sum 
of Diracs. This hypothesis supported by numerical simulations needs however a proof. 



2.75 

2.7 
2.65 

2.6 
2.55 

2.5 
2.45 

2.4 
2.35 



\\u( 




1 / 
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FIGURE 5.7. Evolution of the L 1 -norm of the solution shown in figure 
Left: component u. Right: component w. 
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7. Appendix 



7.1. Derivation of the model. Assume that component v in model (2.5)-(2.7) satisfies 
the steady state equation 



(7.1) 

Solving for v yields 
(7.2) 



= au 2 w + dv — d v. 



a 



d b + d 



u w. 



Substituting ( |7.2[ ) into ( |2.5[ ) and ( |2.7[ ) yields 
(7.3) 



(7.4) 
for cr 
(7.5) 
(7.6) 



u t 



uw 



a- 



a + uw 



d,)u 



-w„ 



7 



d q w — a 1 dbU 2 w + k, 



d b +d 



U; 



After rescaling time, i := d„t yields 



a uw 



d, 



u. 



d g a + uw d g . 

—-w xx - w - a —u w + — 
l d g d g d g 



for x G [0,1], t > 0, 
for x G (0,1), t > 0, 

for x G [0,1], t > 0, 
for x e (0,1). i > 



Defining u(x, t) := J -^j-u(x, t) and w(x, t) 



T £ -w(x,t), we obtain system (2.1)-(2.2): 



(7.7) 
(7.8) 



Wf 



a auw 



d, 



d„ a + cruw d 6 



u 



w — u w + 



ldg ^Jd g d h a 
7.2. Proofs of analytical statements. 



for x G [0,1], i > 0, 
for x G (0,1), i > 0. 



Proof of Theorem 3.1. Since v 



j^u 2 w is the unique root of the right-hand side of 



(2.6), there exists a one-to-one mapping from the set of steady states of (7.3)-(7.4) into 
the set of steady states of (2.5)-(2.7) by (u,w) — > (w, ^^u 2 w,w). 



Since model (2.1)-(2.2) is a linear rescaling resp. linear substitution of (7.3)-(7.4), there 



exists also a one-to-one mapping between the sets of steady states. 

[14] . Theorem 2.6 proves Theorem 3.1 for system ( 2.5[ )-(2.7). Since we found a one-to-one 
mapping between the sets of steady states, statements (ii) and (iii) and existence of the 
steady states in (i) follow from [14], Theorem 2.6. 
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It is left to calculate the exact values of the spatially homogeneous steady states. 
The right-hand-side of (2.1) has two roots: 

(7.9) u = 0, 
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(7.10) 



Ml 



d x 1 

CL\ — d\W 



Substituting 7.9 into the right-hand side of (2.2) and setting it equal to zero leads to 
(7.11) = -tw + Ki, 

defining (u ,w ) = (0,Ki). 



Substituting 7.10 into the right-hand side of (2.2) and setting it equal to zero leads to 

1 



(7.12) 
with roots tZJ_ and w 







-w 



di 



l ai — di w 



+ Ki, 



□ 



To prove 3.2, we use the following lemma from linear algebra, proved in [8], section 
2.1.2: 



Lemma 7.1. Let a real-valued block-matrix 

Ao A Aoo - Dk 2 



(7.13) 



A 



be given with D = diag{d\, d m ),di > 0. 

Let Ai,...,A n denote the eigenvalues of v4 12 and Xi, X m+n the eigenvalues of A. 

Then there exists an injective mapping j : {1, n} — > {1, n+m}, s.t. for all 1 < i < n 

holds 



(7.14) 



lim \jfj\ = Xi, 

k— >oo 



and the real parts of all other eigenvalues of A converge towards — oo as k — > oo. 



Lemma |7.1[ applied to stability of spatially homogeneous steady states of ordinary 
differential equations coupled to reaction-diffusion equations reads: 



Lemma 7.2. Given a system of ordinary /partial- differential equations: 

d 

(7.15) f 

dt 

Let u denote a constant steady state of system ( |7.15 ) and J° denote the Jacobian of the 
ODE subsystem at u: 

d 



-ui = fi(u), l<i<n, 
-Ui = diAui + fi(u), n < i < n + m. 



(7.16) 



J, 



ilu fj(u)\ u=u , 



1 < i < n. 



If J° has a positive eigenvalue X + , the operator resulting from a linearization of (7.15) 
around u has infinitely many positive eigenvalues. 



Proof. The linearization of the right-hand side of (7.15 ) at u — u is of type, written in 
matrix form: 

~J° A 12 
A 21 A 22 -DA ' 



(7.17) 
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and the corresponding eigenvalue problem in matrix form: 
(7.18) 

Assuming <pk being the eigenfunction of the Laplace operator associated to the fcth eigen- 



J°-\ 


A12 




V 




"0" 


A21 


A 2 2 - A - DA 













value, the matrix is of type (7.13). It follows that there exists a sequence of solutions 



(A ( A;), 0^) of the eigenvalue problem (7.18 ) with lim^oo A (fc) = A + , where Re(A + ) > 0. □ 



Now, we can prove Lemma 3.2 



Proof of Lemma 3.2. The Jacobian of the kinetic system of (2.1 )-(2.2) at (u 0l wq) = (0, Ki) 
reads: 



(7.19) 



J 



-di 
-1 



It follows that (0,Ki) is stable solution of (2.1)-(2.2) and its kinetic system. 
The Jacobian of the kinetic system of (2.1)-(2.2) at ( d \ —,,w) reads 



(7.20) 



Since Jn 



J 



(01— jijrfi 



d\ 1 



2di 



V + ((ai-di)tu) 2 ) 



(ai —o!i )d 



is positive, both (u-,W-) and (u+,w+) are unstable solutions of 



(2.1)-(2.2), see Lemma 7.2 To determine stability as steady state of the kinetic system, 



we calculate the determinant and trace of J from (7.20): 

di 



(7.21) |J| 

(7.22) tr(J) 

We note \ J\ ->• - dl(ai ~ dl) 



ax{a\ — di)w 2 
(ai — di)d 



(-(ai- dxfw 2 + dl), 



(ai — <ii)w 



) 2 ). 



«i 



as w — ?• 00. 



The only roots of the determinant |J| are 



(7.23) 



a\ — di 



Since «, ± = f ± J(f )2 - and f > ^ > 0, it follows 



(7.24) 
(7.25) 



|J(u + ,w+)| < 0, 
\J(u-,w_)\ >0, 



what proves instability of (w + ,u> + ), because |J| = A1A2 < 0. 
The stability of (u_,u>_): 

Since J(u_,w^) > 0, (u-^wJ) is unstable if and only if tr( J(w_, wJ)) > 0. 
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tr( J(u-, > is equivalent to 

-(ai - di) - 1 > ' Jl n2 



— (ai - di - 1 > t) -2^> 



^( ai _ dl ) - 1 > (^—^) 2 wl, 



di d 1 f di Y>(—y + 2— /(— ) 2 ( ^ x ) 2 | ^ Kl ^ 2 ' ^ 



ai — di ai — d^ v 2 y 2 V v 2 y K di — d\ J v 2 ai — di ' 



di 



ai — oti zai 2 2 \ 2 a\ — d\ 

This is not satisfied for k 2 > 2 , dl , , . We continue assuming that dl , ^ — (^) 2 > 

± ai(ai— ai) ° ai— ai 2ai \ 2 / 

and define x := ^ and y := dl , . 

2 f ai— ai 

^2 

n' 4 r/ 2 

u l 2 u l 2 ^ 2 2 

—-=y x y > -x y , 

4af ai 

(j/ - ^)x 2 + -^2/ > 0. 
ai 4af 



This is satisfied if and only if 
(7.26) 



dl 

y > — , 

dl 

di 

ai < 



or 



di - 1 

d 4 

%2 < T^yiv-—) 



(7.27) 



d 4 d 2 
Aaf di 



d A 



1 d\ d\ — di{d\ — di) 
Negation yields the result. □ 



Proof of Lemma \3.4\ Adding a multiple of (2.1 ) and (2.2) and integrating over Q leads to 

d f 1 , f ( . u 2 w di , 2 \ 

-u + wdx = / (- u) — w — u w + Ki ax, 



dt J di J \ 1 + uw di 

(7.28) < J (--±u- w + Kl ) dx, 



<— min(di, 1) J ( — u + w^j dx + «i/i(0). 



This leads to 



(7.29) limsup f — / udx + / wdx) < ^ cM^)- 

t^oo \ai J J J mm(d 1 ,l) 
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Additionally, it immediately follows by integrating (2.2) over Q: 



(7.30) 



d 
dt 



J wdx < — J wdx + «;/i(f2). 



From (7.30) follows 



(7.31) 



limsup / wdx < Ki/x(f2). 

t— >oc 



Since w > 0, it follows from (7.28) 



(7.32) 



limsup / udx < — \ r/«i/i(f2). 

t^oo J mm(rfi, 1) 



□ 



Proof of Lemma \3. 5\ The proof is analogues to the proof of Lemma [3T4| without integrat- 
ing over Q. □ 



Proof of Lemma \5.1\ The eigenvector associated to the eigenvalue A + of a 2x2 matrix 

It follows that Vk4>k is the eigenvector of J 



(a-ij + 5 i2 5 j2 D w k 2 ) is v k 
associated to X+(k). 



021 



A+— a.22+D w k 2 



□ 



7.3. Additional figures. In this section, we show numerically obtained solutions which 
were referred to in the previous sections. Additionally, we show for convenience the 



explicit formula for the "perturbation" function p, defined by (5.3)-(5.7): 



4(-l+s-e) 



(s-e)(-2s+2s 2 -e) X 
2{l+2e)x 2 -i(s+t)x+2s 2 +2se-2s 2 e-e 2 

e{~2s+2s 2 -e) ' 
(2s+4s 2 -2s 3 +3t+3st-2s 2 t+e 2 -8x(s+t)+4:X 2 (s+e)) 
(-2s+2s 2 -e)(-l+s+e) : 



(7.33) 



p(x) 



x E [0, s — e), 
x G [s — e, s + e] 
x E (s + e, 1]. 
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Figure 7.1. Left: Finite Fourier Transform f(u) 



Right: Order " o( f'"- 
m_ at a; = 0.250092, x x 



Jg 1 cos(47ra; 2 )e-^ x dx. 



of the growth of the perturbation — e cos(47r;r 2 ) of 

^3 



■r-2 



0.866028 




See fig. 



5.4 



for the 



4 6 8 10 12 14 
t 



Figure 7.2. 



Order u ° ( ? ) - 



-) 



m_ at Si 
solution. 



0.250092 and x 2 



of the growth of perturbation 
See fig. 5.3 



ecos(47rx) of 
for the numerically obtained 
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Figure 7.3. Numerically obtained solution for initial conditions (5.2) with 
e = 0.05, e 1 = 0.1, s = 0.4, 
2.5, d\ = 1.5, Ki = 4 and D w 
w. 



u = 2.215, w = 0.677123 and parameters a x = 
= 5.8541,. Left: component u. Right: component 




Figure 7.4. Numerical solution for initial conditions very close to the stable 
steady state (u ,w ), parameters a x = 2.5,6?! = 1.5, K\ = 4. Left: component 
u. Right: component w. 



7.4. Mesh asymptotic. In this section, we investigate the asymptotic behavior of the 
error due to numerical approximation. Since we do not know the true solution, we inves- 
tigate the asymptotic behavior of the difference of the solution (u, w) and a calculated 
"reference solution" (w re f, w re {). The reference solution is the numerical solution on a much 
finer mesh in time and space. 

First, we show this error for the approximation of the configuration in the introductory 
part for large diffusion coefficient. In that case, only a single spike arises close to the 
position where the initial condition has a maximum. In figure 7J3, the error in L 2 norm 
and the corresponding order of the error reduction under mesh refinement is plotted for 
equidistant mesh. 
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We observe the expected order 0(h 2 ) of error reduction for piecewise linear approxima- 
tion, see e.g. [3]. 



U — Wref £2 \\W — W Te f\\ L 2 




5 10 15 20 25 5 10 15 20 
t t 



Figure 7.5. Upper row: Plot of the evolution of the L 2 -error for a config- 
uration shown in fig. 5.2 and its L 1 norm shown in fig. 5.7 in the sense of 



a reference solution. Lower row: Plot of the evolution of the order of error 
reduction. The reference solution was obtained on a mesh with spatial mesh 
size h = 2~ 13 and temporal mesh size k = 0.01. 
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Figure 7.6. Plot of the evolution of the L 1 -error for a configuration shown in 



fig. 5J2 and its L 1 norm shown in fig. 5/7 in the sense of a reference solution. 
The reference solution was obtained on a mesh with spatial mesh size h = 2~ 13 
and temporal mesh size k = 0.01. 

The same observation holds for the same configuration with smaller diffusion coefficient, 



s.t. growth of more than one spike occurs. The solution is shown in fig. |5.5[ the error in 

fig.O 




FIGURE 7.7. Plot of the evolution of the L 2 -error for a configuration shown 
in fig. 5.5 in the sense of a reference solution. Fig. 5.5 shows the growth of 
multiple spikes due to a smaller diffusion coefficient. The reference solution 
was obtained on a mesh with spatial mesh size h = 2~ 15 and temporal mesh 
size k = 0.00025. 
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It 



Wrcf||L 2 



100 




10 




1 




0.1 


-2- 12 


0.01 


-2- 13 □ ^^<<^^^- 


0.001 


-2~ 14 J^^^j^^ 


0.0001 




le-05 ; 




le-06 < 




le-07 ' 




le-08 


i i i i 



10 



15 



20 



25 




Figure 7.8. Plot of the evolution of the L 2 -error for a configuration shown 
in fig. |5.4| in the sense of a reference solution. Fig. |5.4| shows the growth of 



multiple spikes due multiple maxima of the initial conditions of shape Uq = 
u + cos(27ra; 2 ). The reference solution was obtained on a mesh with spatial 
mesh size h = 2 -16 and temporal mesh size k = 0.00025. 
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s 


■£f=0, max 


•^t=25,max 


0.2 


0.25 


0.2726 


0.4 


0.417 


0.43237 


0.5 


0.5 


0.5 


0.7 


0.66 


0.645 


0.85 


0.792 


0.77 



Table 1. Position x 
with e = 0.1, ei 



max of the arising spike (t = 25) for initial conditions 
0.05, D w = 6 and maximum at x t= o, m ax- The shape 



of solutions are as in Figure |5.2[ differing qualitatively only in the position 
of spike/sink. We observe that a spike grows close to the position of the 
maximum of the initial conditions. 





spikes 


D wA = 5.8541 

I n 
l n 

l n 

25 

l n 

3fi f jl 


1 

2 
3 
3 
4 
4 



Table 2. Number of spikes arising for different diffusion coefficients Di, 
initial conditions (5.2) with s = 0.4, ei = 0.05, e = 0.1. 
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